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ABSTRACT 



We calculate the distribution of neutral Hydrogen within 750 proper h ^kpc of a 
quasar, Lboi = 1-62 x IO^^Lq « ^Edd, powered by accretion onto a super massive black 
hole, Mbh = 4.47x IO^Mq, at z = 3. Our numerical model includes cosmological initial 
conditions, gas dynamics, star formation, supernovae feedback, and the self consistent 
growth and thermal feedback of black holes calculated using GADGET as well as a 
detailed post-processing ray tracing treatment of the non-uniform ionizing radiation 
field calculated using SPHRAY . Our radiative transfer scheme naturally accounts for 
the self shielding of optically thick systems near the luminous central source. We show 
that the correct treatment of self shielding introduces a flattening feature into the 
neutral column density distribution around Log iVjji = 20 and that regions with the 
lowest neutral fractions are not necessarily those with the highest density gas. 

For comparison with our numerical work, we solve a Ricatti equation which deter- 
mines the equilibrium Hydrogen ionization fractions in the presence of a radiation field 
that falls off as one over r squared with regions above a given gas density threshold 
completely shielded from ionizing radiation. We demonstrate that these simple semi 
analytic models cannot reproduce the neutral Hydrogen field calculated using SPHRAY. 
We conclude by comparing our models of this single proximity zone to observations 
by Hennawi and Prochaska of the absorption spectra of background quasars which 
are coincident on the sky with foreground quasars in their Quasars Probing Quasars 
(QPQ) series of papers. Compared to the QPQ sample, we find a factor of 3 fewer 
optically thick (Log A^hi ^ 17.2) systems around our quasar, however the dark matter 
halo that hosts our simulated quasar, Mnaio = 5.25 x IO^^Mq, is less massive than 
the typical QPQ host halo by a factor of four. Allowing for a linear scaling between 
halo mass, baryonic overdensity and number of absorbers, we estimate the typical host 
halo mass in the QPQ sample as 1.92 x IO^^Mq. 

Key words: astrophysics, theory, numerical, simulation, SPH, ray tracing, simula- 
tion, radiative transfer, quasar, AGN, black hole 



1 INTRODUCTION 

After the reionization of the Universe by the first luminous 
objects, the majority of neutral Hydrogen resides in gravi- 
tationally collapsed objects, specifically Damped Lyman-Q 
systems (DLAs). These systems are observed via absorp- 
tion lines in the spectra of distant quasars and are histor- 
ically identified as those absorbers with HI column densi- 
ties A'hi > 2 X lO^^cm"^ with lower column density sys- 
tems 10^^'^ < iVni < 2 X 10^°cm~^ being labeled Lyman 
Limit Systems (LLSs, see lWolfe et al.l i2005. for a review). In 
contrast, the systems that gi ve rise to the Lyman-g forest 
have A'hi < lO^'^'^cm"'^ (see lRaucb|[l99i ; [Weinberg et al.l 
|2003| . for reviews). This historical column density thresh- 
old (A'hi ~ lO^^'^cm"^) serves to divide systems into those 



that are predominantly neutral (the DLAs and LLSs) and 
those that reside in the mostly ionized intergalactic medium 
(IGM). The DLAs and LLSs remain mostly neutral in the 
presence of an ionizing ultraviolet (UV) background and lo- 
cal point sources through self shielding. Studying them in 
absorption opens a window into the post reionization pop- 
ulation of cold, dense, neutral gas and provides a survey 
technique with a bias complimentary to that of emission 
surveys. 

To study these systems numerically requires treat- 
ing the UV ionizing background at some level. Many 
cosmological simulation s h ave followed in th e steps of 
iHaardt fc Madaul (|l99i) and lKatz et"al] (llQeeal l by consid- 
ering a spatially uniform, time variable, background with 
a spectral shape characteristic of quasar and stellar radia- 
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tion that has been reprocessed by the IGM. The background 
plays a role in determining the cooling function of cosmolog- 
ical gas and so its inclusion is also necessary at some level 
for a realistic description of galaxy and star formation. A 
uniform background is a good first approximation if one is 
interested in radiative cooling, however it ignores completely 
the self shielding that defines the DLAs and LLSs. 

The most straight forward way to account for this 
shielding is to reduce the UV field to zero i n regions above 
a give n gas density threshold as was done in iHaehnelt et alj 
l|l998t ). A more detailed treatment based on the solution for 
plane parallel radiation in c ident o n a constant density slab 
is described in iKatz et al.' (■ 1996b') and used in the work of 
[Gardner etlo] l|l997al . |2001„ 1997b l . These corrections are 
not based on transferring radiation through the simulation 
volume, but rather on applying the plane parallel solution 
on a pixel by pixel basis to HI column density maps in post 
processing. 

Other appr oaches to self shielding include that of 
ICen et al] (|2003l ) who include attenuat ion of the UV back - 
ground on a cell by ce ll basis akin to iKatz et al.l (Il996d ). 
iRazoumov et al] (|2006l ). who include a treatment of the 
transfer of ionizing radiation in post processing using the 
Fully Threaded Transport Engin e (FTTE) descr i bed i n 
IRazoumov fc Cardalll l|2005l '). And iNagamine et al.l l|2007l ). 
who use a multi phase gas model to treat star formation 
and assume the cold dense phase to be fully neutral. 

Recently, a series of nested galaxy formation simu l ations 
at different resolutions was used by IPont zen et alj (|2008l ) 
to study DLAs over a broad range of mass scales. They 
include a simplified radiative transfer scheme to account for 
self shielding which lies somewhere between the full radiative 
transfer modeling and the pixel by pixel corrections. 

The works mentioned above were aimed at studying sys- 
tems where the radiation field is n ot dominated by a local 
point source. iMiralda-Escudel (|2005l ) applied simple analytic 
arguments based on the conservation of surface brightness 
to argue that on average the effect from local sources is neg- 
ligible compared to that of the background for systems with 
optical dept hs below that of Lyman Limit Systems. On the 
other hand, ISchavd (|2006l ') use analytic arguments to show 
the local radiation field is likely to be important for denser 
systems that tend to cluster around the large scale overden- 
sities that host these sources. 

In this work, we examine the balance between the UV 
background, local sources, and self shielding by combining, 
for the first time, a hydrodynamic simulation that tracks 
the formation and accretion history of black holes with a 
detailed ray tracing treatment of the non-uniform UV ra- 
diation field. Shielding is most important in the presence 
of dense gas and strong UV fields, both of which occur 
near active galactic nuclei (AGN) and quasars. In fact, the 
transition between the background UV field and the local 
AGN/quasar UV field serves to define the proximity region 
of these objects. With the availability of large quasar cata- 
logues and high resolution spectroscopy, it has become pos- 
sible to study absorbers proximate to a foreground quasar in 
the spectrum of a coincident background quasar. Work such 
as this has been carried out in a series of papers entitled 
"Q uasars Probing Quasars", the first of which is authored 
bv lHennawi et al] l|2006l ) (HP06 in the rest of this paper). 
We compare our numerical models to observations from this 



body of work a nd to theoretical calculati o ns of the HI col- 
umn density bv lZheng fc Miralda-Escudel l|2002l ). 

The format of this paper is as follows. In §2 we de- 
scribe the GADGET simulation that determines our tem- 
perature and density fields, in §3 we review the ray tracing 
code SPHRAY used to calculate the transfer of ionizing ra- 
diation through this density field, in §4 we describe how we 
model our sources of ionizing radiation, in §5 we describe our 
semi analytic model, in §6 we describe our ray tracing re- 
sults and compare them to the semi analytic model, in §7 we 
compare our theoretical results with those of iHennawi et al] 
(2006), and IZheng fc Miralda-Escudel l|2002D and in §8 we 
conclude and discuss further work. 



2 COSMOLOGICAL HYDRODYNAMIC 
SIMULATION 

The proximity zone we study in this work was cut from a 
cosmological simulation performed using a modified version 
of the publicly available S moothed Partic le Hydrodynam- 
ics (SPH) code GA DGET (jSpringell [20051'). The interested 
reader is directed to iDi Matteo et aL ( 20081 ) for a full de- 
scription of the simulation, however we will briefiy describe 
it here. 

We have adopted cosmological p arameters consisten t 
with the WMAP first year results jSpergel et alj |2003| ). 
specifically {Ho,Q,m,^A,^h,^s) = (70,0.3,0.7,0.04,0.9), 
where Ho is the Hubble parameter today, f2m, f^A, and fib 
are the total matter, dark energy, and baryonic matter den- 
sity parameters, and as is the mass variance in spheres 
of 8/i~^Mpc. The simulation volume is a periodic cube 
of side length 33.75 comoving h~^Mpc. The matter dis- 
tribution is sampled initially □ by 486^ dark matter and 
486"^ baryonic particles, resulting in a mass per particle of 
2.41 X 10^ /i~^MQ for the dark matter and 3.72 x 10** fe'^M© 
for the baryons. The simulation was evolved from redshift 
= 99 to 2 = 1, and we centered our cutout on the most 
luminous accreting black hole from a snapshot at z = 3 near 
the peak of quasar activity. 

The main physics elements o f the simula tion include 
a Tree-PM solver for gravity (e.g. lBaglall2002l ). an entropy 
conserving implementation of SPH for hydrodynamics from 
ISpringel fc Hernguistj l|2002l) . a model for star formation and 
feedback by supernovae described in ISpringel fc Hernguistl 
(l2003h . optically thin radiative cooling in the presence of a 
spatially uniform UV background as in iKatz et all (|l996al ). 
and a model for the accretion of gas onto bl ack holes and the 
associated thermal feedback de scribed in iDi Matteo et al] 
(|2005h and ISpringel et al.l (|2005| ). 

This GADGET simulation provides the base density 
and temperature fields for our work. In Figure [1] we show 
their distribution. The red grouping in the lower left corner 
of each panel represents the underdense to mildly overdense 
intergalactic medium (IGM) responsible for the Lyman a 
forest. Its slope determines the effective equation of state 
of the IGM. The green/yellow grouping above that is the 



^ During the course of the simulation, baryonic particles can lose 
mass by spawning stellar particles and can disappear all together 
through accretion onto black hole sink particles. 
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Figure 1. Density vs Temperature distribution. These quantities 
are fixed across all radiative models. Shown in the top (bottom) 
panel is the mass (volume) weighted probability density distribu- 
tion in Log A — T space where A is the gas density in units of 
the baryon critical density and T is the temperature in Kelvin 
[i.e. an integral over the plane shown will give unity]. The color 
bar is logarithmic such that a unit difference between two pixels 
indicates a factor of ten difference in mass (volume). 

warm/hot intergalactic medium. The check mark shaped re- 
gion above A « 1000 in the mass weighted panel represents 
collapsed objects where the bulk of neutral Hydrogen re- 
sides. 



3 RADIATIVE TRANSFER SIMULATIONS 

The thermal state of the gas in the hydrodynamic GADGET 
simulation described above is set by photo heating from a 
uniform UV background, shock heating, cosmic expansion, 
and direct energy injection due to feedback from supernovae 
and accreting black holes. The ionization state is set by as- 
suming photo ionization equilibrium with the uniform UV 
background, and taking account of the temperature depen- 
dant coUisional ionization and recombination rates. We wish 



to improve the treatment of the UV field while preserving 
the local thermal feedback effects from supernovae and black 
holes. To do this, we select the black hole with the largest 
accretion rate from a GADGET snapshot at z = 3, cut 
out a 6/i^^Mpc (comoving) region centered on it, and fix 
the density and temperature fields. We then calculate the 
equilibrium ionization state of the gas in the presence of a 
non uniform ionizing radiation field. To do this, we use two 
methods. The first is a semi analytic model which includes 
a self shielding gas density threshold above which the pho- 
toionization rate is set to zero. The second is a detailed ray 
tracing calculation of the non uniform UV radiation field 
which naturally accounts for self shielding. This ray tracing 
is accomplished with SPHRAY . 

A deta i led d escription of the code can be found in 
lAltav et al.l (|2008t ). however we will outline the important 
features here. SPHRAY is a long characteristics, monte carlo, 
ray tracer for radiative transfer post processing. It works by 
transporting a large number of photon packets through the 
smoothing kernels of an SPH density field. The origin, direc- 
tion, and frequency of each packet is sampled from user sup- 
plied probability distribution functions. It does not require 
that the SPH density field be smoothed and so can preserve 
its Lagrangian and adaptive nature. For this project, it was 
necessary to update the ionization solver to handle the very 
small neutral fractions (on the order of 10~*). 

For all of the radiative transfer (RT) calculations in this 
paper, the Hydrogen ionization fractions and temperatures 
of all particles are initialized to the values in the GADGET 
snapshot. The density field is then exposed to one of sev- 
eral non uniform radiation fields (see §4) until it reaches an 
equilibrium state. In the following sections, we make com- 
parisons between the neutral density fields produced using 
SPHRAY and the semi analytic model. 

We are justified in calculating the equilibrium neutral 
fractions by the fact that at z = 3 we are not calculating 
the expansion of an ionization front but the response of the 
gas to an inhomogeneous radiation field that has been in 
place long enough for the gas to respond. The light crossing 
time from the central quasar to the corner of the simulation 
volume is « 6 Myr. This means that light travel times will 
have to be taken into account if we want to consider variable 
luminosities on time scales shorter than this. For the current 
work we assume the central point source has had a steady 
luminosity long enough for the influence of this source to 
propagate to the whole cut out and for the gas in the cut 
out region to respond to this radiation fleld. 

We have used the case A recombination rates as most 
of the photons emitted by the recombining gas will leave 
the highly ionized computational volume. This is not the 
correct approximation for the self shielded gas, however the 
very small ionization fractions in these collapsed objects im- 
ply that only a small fraction of the gas there is actually 
recombining. 



4 SOURCES OF IONIZING RADIATION 

In this section we describe our various source models. We 
decompose the radiation fleld into components due to the 
UV background which comes into the simulation volume 
from the boundaries, the central quasar, the smaller satel- 
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Table 1 . Summary of the Sources of Ionizing Radiation 
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Figure 2. Luminosity of the accreting black holes contained 
within our cutout volume vs their mass. The bottom panel shows 
the number of Hydrogen ionizing photons emitted per second as- 
suming a spectral slope of auv = 1-76. This plot takes into ac- 
count bolometric corrections. Also shown in the bottom panel 
(dotted line) is the contribution from the background UV field. 
The top panel shows the bolometric luminosity of each black hole 
and lines indicating the Eddington luminosity, I/Bddi a^nd one 
hundredth LEdd- The central quasar in our cutout region (the 
most luminous source) is the filled circle in both panels. 
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lite accreting black holes, and the star forming galaxies. We 
have summarized the fiducial luminosities in Table [l] Note 
that 98.4% of the ionizing photon flux comes from the back- 
ground and the central black hole. There is a 1.60% contri- 
bution from the minor black holes and a 0.03% contribution 
from the young stellar clusters. For this reason, we have cho- 
sen to include only the background and the central source in 
our models. Below we describe how we calculate the photon 
flux from these sources. 

We assume that, above the HI ionizing threshold, the 
spectra of all our source populations can be modeled with 
power laws of the form. 



J = Jm — , (1) 

where J is a specific intensity with units of ergs cm"'^ s~^ 
Hz^^ sr-\ 
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Figure 3. Projection of the gas density and source positions in 
the cutout region onto the x-y plane. The image was made by col- 
lapsing the 750^/i~'^kpc^ volume along the z-axis. Shown in black 
and white is Log A where A is the average gas density in each 
column in units of the critical gas density. The red circles repre- 
sent accreting black holes. The size of each symbol is proportional 
to the logarithm of the flux. The circle in the bottom left corner 
represents the magnitude of the background flux coming in from 
outside the cut out region. The radiation field is dominated by 
the central black hole and the background flux. 



4.1 Accreting Black Holes 

The GADGET simulation tracks the formation and ac- 
cretion history of black holes. The bolometric luminosity 
of each black hole is related to its accretion rate through 
an efficiency parameter r\ = Lh/M-BYi(?- We take as our 
fiducial value = 0.1 which is the mean value for radia- 
tively efficient accretion on to a Schwarzschild black hole 
(jShakura fc Svunvaev|[l973h . In S7 we also consider models 
with smaller and larger radiative efficiencies. In the GAD- 
GET simulation a small fraction of this energy is thermally 
coupled to the nearby gas E^ccd = ^fLb- A value of e/ = 0.05 
brings the simulation into agreement with obser vations of 
the Mbh - o relationship l|Di Matteo et al-lliooH ). 

This population of sources consists of a dominant cen- 
tral black hole and tens of black holes with smaller accretion 
rates surrounding it (see Figures [2] and |3]). For each of these 
black holes we calcul ate a B band luminos ity using the bolo- 
metric corrections of lMarconi et al.l (|2004l ). To translate this 
to a luminosity at the Hydrogen ionizing threshold we use a 
broken power law spectrum with spectral index aop = 0.44 
in the optical part of the spectrum and auv = 1.76 in the 
UV part with the break at 1200 Angstroms. Given the lu- 
minosity at the ionizing threshold, we integrate the power 
law from 1 to 36 Rydbergs to calculate the number of ion- 
izing photons per second. The choice of cutoff frequency is 
dictated by a deviation from the UV power law above ~ 36 
Rydbergs in the quasar template spectra and the fact that 
high energy X-ray photons above 36 Rydbergs have mean 
free paths much longer then our box size. 

With a box of this size, it is not possible to examine the 
full quas ar luminosity function due to the lack of rare objects 
however, iDeGraf et aP (|2009l ) examine its faint end using 
the base hydrodynamic simulation discussed above and two 
similar simulations. 
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4.2 UV Background 



5 IONIZING RADIATION MODELS 



The ionizing background is a combination of the ionizing ra- 
diation produced by quasars and star forming galaxies that 
reside outside our simulation volume. Estimating the magni- 
tude of this radiation field observationally is challenging. In 
general, there have been two approaches. One involves inte- 
grating observed luminosity functions and another involves 
inferring the background from the HI photoionization rate 
calculated using quasar absorption lines. The latter method 
has the advantage that the absorption lines probe the local 
z = 3 gas. 

The fundamental observable we've used to calibrate the 
background flux used in our simulation is th e Lym an a for- 
est measurements of iFaucher-Giguere et al.l (|2008l ) (FG08) 
who calculate a mean Hydrogen photoionization rate Fhi = 
(0.59 ± 0.07) X 10"^^ s"^ at z=3. The Hydrogen photoion- 
ization rate for a uniform background with spectral index a 
is given by, 



Thi = 47r / —am Jhi 



(2) 



where h is Planck's constant, cthi is the photoionization cross 
section of Hydrogen at the Lyman limit, and Jhi is the angle 
averaged specific intensity. Integrating and solving for Jhi 
we get. 



Jhi = 



Tmhia + 3) 

47r(THI 



(3) 



Using the slope we adopt for our background spectrum 
Q = 1.00 gives a numerical value of Jhi = 2.35 x 
10"^^ ergsHz-^s"^ cm and implies a photon number den- 
sity, 



47r 



hv 



-dv 



(4) 



For a uniform UV background along with the spectral 
index a completely characterizes the radiation field. 

The emissivity e, of the sources is connected to the 
background intensity Jhi through the mean free path /hi 
of the photons. In the Local Source approximation, sources 
separated from a given point by more than one mean 
free path have their con tribution highly attenuated (see 
ISchirber fc BuUockl l20o3) . In this regime the emissivity is 
expressed simply as. 



e = 47r 



Jhi 



(5) 



This m ean free path is calculated in iFaucher-Giguere et al] 
l|2008h as 



Im = 85 



l + z 



proper Mpc 



(6) 



At 2 = 3, approximately half of the ionizing photons in the 
background originate in star forming regions. This emissivity 
and the relatively small volume of our simulation compared 
to Zhi imply that less than 1 in 1000 ionizing photons we 
consider originate in star forming regions within the simu- 
lation. 



The ionization state of Hydrogen is a balance between 
photo/collisional ionizations and recombinations. In partic- 
ular, the Hydrogen ionization fraction x = nnii/wH, evolves 
according to the equation 



dx 
H 



[Fhi + 7(r) n4x)]x - a(T) n,{x) x, 



(7) 



where Fhi is the photoionization rate, 7 is the coUisional 
ionization rate, a is the recombination rate, and rio is the 
electron number density. If we decompose rie into a part 
due to Hydrogen ionizations and a part due to Helium and 
metal ionizations rie = (x -\- y)nH we can rewrite the above 
equation as. 



dx 
It 
Rit) 

Q{t) 
Pit) 



= Rx^ + Qx + P, 

= (7 + a)nH 

= Fhi — 7nH + (7 + a)nHy 

= -(Fhi+7^h2/) 



(8) 



This is a Riccati equation and the equilibrium solution 
is given by the positive root of the quadratic formula. In 
a pure Hydrogen and Helium gas, the variable y is bound 
between y = for completely neutral Helium and y = (1 — 
X)/2X for fully ionized Helium, where X is the Hydrogen 
mass fraction of the gas. 



5.1 Semi Analytic Model 

For comparison with our ray tracing results, we calculate 
the neutral Hydrogen field using several simplifying assump- 
tions. Given a density field nn, a temperature field T, an 
electron abundance from Helium and metals y, and a pho- 
toionization rate Fhi, we can analytically calculate the equi- 
librium Hydrogen ionization fractions using equation (IS}. 
We take the first three quantities from the base hydrody- 
namic run and supply the photoionization rate by hand. 

The simplest model contains only the unshielded back- 
ground UV field. This amounts to setting Fhi = ^m^'^ ~ 
0.59 X lO"^'^ for each particle. To include the central 
source we add to this background an inverse squared contri- 
bution from the quasar. 



nQSO 



1 du (ji_y\ (j^ 

47rr2 I hv \ imi J \ vm 



' ■'HI 

^hichi 



(9) 



47rr2/i(auv + 3) 

where Lhi is the luminosity density of the quasar at one 
Rydberg in units of ergs s~^ Hz~^ and quv is the spectral 
slope in the UV part of the quasar spectrum. For these mod- 
els the photoionization rate Fhi ~ r^f"'' + F^^'^(r). Finally, 
we include a set of models in which Fhi = F^f"'' -I- F§^'^(r) 
for particles below a given gas density threshold and zero 
above it. These models are labeled by the density threshold 
in units of the critical gas density on a logarithmic scale (e.g. 
"A = 4.0" has gas with density greater than ten thousand 
times the mean shielded). 
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Figure 4. Mass weighted probability distribution of Fhi for sev- 
eral background only ray tracing models. The arrow indicates 
the observed value at 2: = 3 from FG08. The width of the ar- 
row head indicates the l-tr error bars on their measurement. For 
a completely uniform background this distribution would be a 
delta function at the arrows position. The primary peaks are due 
to the mildly over dense IGM while the smaller secondary peaks 
and tails towards small Fjji are due to self shielding. 

5.2 Ray Tracing Models 

For these models we calculate the transfer of ionizing radi- 
ation emitted from the central quasar inside our computa- 
tional volume as well as a background flux sent in through 
the boundaries. These models naturally account for the self 
shielding of dense particles. The flux that constitutes the 
background is calibrated by running several ray tracing mod- 
els without the central quasar contribution. The flux that 
most nearly reproduces the Fhi observed by FG08 is set 
as the fiducial background strength (see Figure Q. This 
flux is used in all the models that include the contribution 
from the central quasar. The luminosity of the black hole 
is determined by the accretion rate from the base hydrody- 
namic simulation. We have introduced an arbitrary scaling 
factor into the central quasar luminosity to model the ef- 
fect of changing the efficiency with which accreted matter is 
converted into ionizing photons. These models are labeled 
"GADGET -I- SPHRAY" . 



6 RESULTS 

In this section we discuss the radiative transfer results and 
make our comparisons between the ray tracing and semi 
analytic models. 

6.1 Neutral Fractions 

First we will take the global view and examine the distribu- 
tion of neutral fractions for all the particles in our simula- 
tion volume. In Figure [S] we show this distribution for sev- 



Figure 5. Mass weighted probability distribution of xhi for 
several background only models. The dashed green line indicates 
the distribution taken directly from the GADGET simulation. 
The solid black line indicates the distribution after the GAD- 
GET data has been post processed with SPHRAY using the 
fiducial background strength. The dashed blue line is a semi an- 
alytic model in which the photoionization rates were set to zero 
everywhere and only coUisional ionizations were considered. The 
red dash-dot line is a semi analytic model with the shielding den- 
sity threshold set to infinity (i.e. zero shielding) while the orange 
dash-dot-dot-dot line is identical except that the shielding thresh- 
old has been set to ten thousand times the mean density. The 
total amount of neutral Hydrogen in each model is indicated in 
the legend. 



eral models that include only a UV background. The triple 
peaked structure of this distribution is common to all mod- 
els, the left peak being due mostly to coUisional ionizations, 
the central peak due to photoionizations in the IGM and the 
right peak formed by particles self shielding from the ion- 
izing radiation. The vast majority of the neutral material 
resides in the self shielding peak and we have indicated in 
the Figure the total amount of neutral Hydrogen for each 
model. We present these background only results to demon- 
strate the consistency between the original GADGET sim- 
ulation, the radiative transfer post processing results, and 
the semi analytic models away from the self shielding re- 
gion. This is the reason why a detailed treatment of self 
shielding is not necessary for studies of the Lyman-a forest, 
but is necessary for studies of denser systems. 

In Figure |6] we show the xm distribution for several 
radiation models that include not only a UV background 
but the central quasar as well. The results directly from the 
GADGET simulation and those that include only coUisional 
ionizations are reproduced from Figure [5] for reference. The 
semi analytic model in which particles with a total gas den- 
sity above Log A = 4.0 are completely shielded reproduces 
the total amount of neutral Hydrogen present in the fiducial 
ray tracing model, however the self shielding peak is artifi- 
cially concentrated at very high neutral fractions. We will 
see how this difference affects the radial HI profiles in the 



HI in QSO Proximity Zone 7 



0.35 



0.30 



0.25 



0.20 



0.15 



0.10 



0.05 



0.00 





Godget 

Godget + Sphroy 

Collisionol Only 

Anolytic ^ = CO 

Ar^lytic A = 4.0 



—1 I I L 



1.5 



2.0 

Log R [kpc/h] 



2.5 



Figure 6. Mass weighted probability distribution of xhi for sev- 
eral background plus quasar models. The solid black line indicates 
the distribution after the GADGET data has been post pro- 
cessed with SPHRAY using the fiducial background and quasar 
luminosities. The red dash-dot line is a semi analytic model with 
the shielding density threshold set to infinity (i.e. zero shielding) 
while the orange dash-dot-dot-dot line is identical except that the 
shielding threshold has been set to ten thousand times the mean 
density. The dashed green line indicates the distribution taken 
directly from the GADGET simulation which included only a 
uniform UV background. The dashed blue line is a semi analytic 
model in which photoionization rates were set to zero and only 
coUisional ionizations were considered. These last two models do 
not include a central quasar component (or any photoionizations 
in the latter case) and are shown only for reference. The total 
amount of neutral Hydrogen in each model is indicated in the 
legend. 



Figure 7. Cumulative HI distribution. Shown is the neutral Hy- 
drogen mass [MfTj/h] within a given radius R [proper kpc/h] for 
several background plus quasar models. The solid black line indi- 
cates the distribution after the GADGET data has been post pro- 
cessed with SPHRAY using the fiducial background and quasar 
luminosities. The red dash-dot line is a semi analytic model with 
the shielding density threshold set to infinity (i.e. zero shielding) 
while the orange dash-dot-dot-dot line is identical except that the 
shielding threshold has been set to ten thousand times the mean 
density. The dashed green line indicates the distribution taken 
directly from the GADGET simulation which included only a 
uniform UV background. The dashed blue line is a semi analytic 
model in which photoionization rates were set to zero and only 
coUisional ionizations were considered. These last two models do 
not include a central quasar component (or any photoionizations 
in the latter case) and are shown only for reference. We have also 
added several lines that go as for reference. 



next section, the A vs. xm distribution in § 6.3 and the HI 
column densities in § 6.4. The completely unshielded semi 
analytic model produces no particles with neutral fractions 
above xm ~ 0.02 and contains only 2 % of the neutral Hy- 
drogen present in the fiducial ray tracing model. 



6.2 Radial Distributions 

Here we show how the neutral fraction PDFs translate into 
radial neutral mass distributions. In Figure [T] we plot the 
cumulative HI mass within a given radius for the same mod- 
els as in Figure |6] As was mentioned in the previous section, 
the semi analytic model in which particles with a gas density 
Log A > 4.0 are shielded yields the same total neutral mass 
as the fiducial ray tracing model, but the radial profiles are 
significantly difi'erent. This means that self shielding efi'ects 
are important not only for absorption line measurements but 
also for HI emission observations attempting to connect the 
HI mass with the underlying total density field. 



6.3 Density vs Neutral Fraction 

Here we examine the relation between neutral fraction and 
total gas density. We show this distribution for three semi 
analytic models and our ray tracing model in Figure H) 
Naively, one might associate the densest gas with the most 
neutral gas, however this is not the case due to the high 
temperatures produced by shock heating and feedback which 
cause coUisional ionizations. Instead, most of the neutral ma- 
terial is in a peak at approximately 10,000 times the mean 
density where shielding can occur and temperatures remain 
in the tens of thousands of degrees. 

In the third panel we show the results for a semi analytic 
model with a shielding threshold at Log A = 2.8. We include 
this panel not because this is a realistic density threshold, 
but rather to show how the distribution above the threshold 
compares to the ray tracing results in the last panel. 

6.4 Neutral Hydrogen Column Densities 

We have calculated the distribution of neutral Hydrogen col- 
umn densities through the 1500 kpc/h (physical) simulation 
volume. To do this, we collapsed all the gas along the z-axis 
and then calculated the column density on each of 2048 x 
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Figure 8. Density vs Neutral Fraction distribution for several radiative models. All of them include the fiducial strength background 
and central quasar. The panel on the left is a semi analytic model with no shielding, the next two are semi analytic models with shielding 
thresholds at A = 4.0 and 2.8. The last panel is the fiducial ray tracing model. Shown is the neutral mass probability density distribution 
in Log A — zjji space where A is the gas density in units of the baryon critical density and xjji is the neutral fraction [i.e. an integral 
over the plane shown will give unity] . The color bar is logarithmic such that a unit difference between two pixels indicates a factor of ton 
difference in neutral mass. 




Log(N„|[cm^]) 

Figure 9. Distribution of Neutral Hydrogen Column densities 
obtained by collapsing the 1500 kpc/h simulation volume along 
the z axis. The top panel includes all lines of sight with an impact 
parameter b < 750 kpc/h, while the bottom panel is restricted to 
those lines of sight with an impact parameter 6 < 150 kpc/h. The 
GADGET model does not include the UV field from the central 
quasar. We have indicated the maximum y-value in the upper 
panel with a horizontal line in the lower panel. 



2048 pixels on this collapsed plane interpolating from the 
SPH smoothing kernels. This leads to 0.732 pc/h between 
hypothetical lines of sight. 

In Figure [9] we show the probability distribution func- 
tion of column densities for all the lines of sight with impact 
parameters b < 750 kpc/h (top panel) and b < 150 kpc/h 



(bottom panel). The peak of this distribution is off the plot 
to the left as most of the area on the collapsed plane is in 
the IGM, however in this work we are concerned with the 
optically thick systems above Log Nhi ~ 17.2. 

We examine this distribution in four models. Three of 
these include the central quasar and background UV fields 
and one, the direct GADGET output, only includes a back- 
ground UV field. This lack of a central source in the latter 
model causes the PDFs in the top and bottom panels to have 
the same shape as the restriction to small impact parame- 
ters does not select a region of space with an enhanced UV 
field. It is interesting to note that the A — 4.0 semi analytic 
model shows this same radial independence. Examining Fig- 
ures [7] and 1141 it is clear that the semi analytic model with 
zero shielding and the ray tracing GADGET -I- SPHRAY 
model have cleared almost all the HI from the central 100 
kpc/h. In the A = 4.0 model, this HI is largely filled back 
in which removes the radial dependence. This highlights the 
main problem of using a strict density threshold criteria for 
semi analytic self shielding in the presence of point sources. 
In this case it is clear that this choice of threshold produces 
too much shielding close to the source and too little shielding 
at larger radii. 

In the top panel of Figure [51 each model demon- 
strates, albeit to varying degrees, a power law at low col- 
umn densities followed by a flattening and then a cut off 
at higher column densitie s . This behavior was examined by 
I Zheng fc Miralda-Escudel (|2002l ) (see Figure 2 in that work 
for their Nhi distributions) using a model in which a spher- 
ical and isothermal gaseous halo is illuminated by a cosmic 
background UV field and ascribed the flattening feature to 
self shielding. Our ray tracing model, with the most realistic 
treatment of self shielding, produces the largest flat section, 
however the semi analytic A = oo model and the direct 
GADGET output also contain slight flattening features while 
neither model includes any self shielding. In addition, this 
flattening feature is washed out near the quasar for the semi 
analytic A = oo and the ray tracing model (bottom panel). 
This means the intrinsic density field is such that it pro- 
duces a slight flattening feature and in the presence of mild 
UV fields this feature is enhanced by self shielding. How- 
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ever, once the UV field becomes too strong, the feature is 
completely washed away. 

A simple statistic to derive from these data is the cov- 
ering fraction of pixels above a given column density thresh- 
old /i^^' defined as the number of pixels above A''hi divided 
by the total number of pixels. In Table 2, we report these 
covering fractions fc for lines of sight with Log Nhi above 
17.2 and 19.0 for the same impact parameter cut ofTs as 
above. We also inclu de the covering fractions reported in 
iHennawi etlo] (|2006l '). We aUow for various central quasar 
luminosities in our ray tracing models and several shielding 
thresholds in our semi analytic models. These numbers do 
not take into account the fact that the observational sample 
of HP06 does not uniformly cover impact parameter space, 
but we address this issue in §7. One can see from the table 
that shielded semi analytic models that manage to repro- 
duce the covering fraction from the ray tracing models do 
not reproduce the total HI mass in the cylinder. One can 
also see that both of our models fall short of the covering 
fractions implied by HP06. There are several explanations 
for this and they will be discussed fully in §7. 



6.5 Photoionization Rate 

We examine the radial and angular dependence of the pho- 
toionization rate from our fiducial GADGET + SPHRAY 
model in Figures \W\ and [11] The PDFs in Figure [10] serve 
to illustrate two things. The equal spacing between peaks 
demonstrates the general inverse square character of Fhi 
and the width of the distributions on each shell show the 
effects of partial shielding from the ray tracing run. This is 
far from a binary "on" or "off" distribution. We note that 
several authors have utilized more detailed semi analytic 
prescriptions than ours for self shielding, however without a 
detailed radiative transfer model to compare with, the suc- 
cess of these models at reproducing the true shielding is in 
question. We have indicated the photoionization rate which 
an unshielded particle would feel at each radii using vertical 
lines. 

To visualize the angular dependence of Fhi on these 
shells in the fiducial ray tracing run, we use the HealPix 
algorithm to create MoUweide projections of this quantity 
at the same radii as in Figure [TD] The self shielded particles 
show up as dark spots and some shadowing is revealed as a 
persistent dark spot in the lower left of each diagram. 

A radiation field that is restricted to limited solid an- 
gles would help explain the non-detection of the transverse 
proximity effect (e.g.[Croft 2004). Because large dark matter 
halos which host quasars typically lie at the intersection of 
several filaments, one might expect that isotropic emission 
could be reduced in some solid angles by filtering through 
the surrounding density field. We have shown here that this 
is not the case and that if the radiation from quasars is re- 
stricted to certain solid angles that obscuration by optically 
thick material near the black hole (below the resolution of 
our simulation) is a more likely explanation. 



7 COMPARISON TO OBSERVATIONS 

In the fir st of a series of paper s entitled "Quasars Probing 
Quasars" iHennawi et al.l (|2006l ) have compiled a catalogue 
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Figure 10. Probability distribution function of FHI on a series 
of spherical shells. The radii are R = 12.5, 25, 50, 100, 200, and 
400 kpc/h from the central black hole. The equal spacing of the 
peaks indicates a general behavior and the width indicates 
shadowing and shielding. We have indicated the photoionization 
rate which an unshielded particle would feel at each radii using 
vertical lines. 



of quasars with small angular separations on the sky. This 
sample can be used to study the environment of the fore- 
ground quasars using the absorption lines in the spectra 
of the background quasars. Specifically, they measure the 
number of sight lines that intersect absorbers with column 
densities above Log A^hi ~ 17.2 and Log Ahi ~ 19 corre- 
sponding to Lyman Limit systems and Super Lyman Limit 
Systems. 

We now make a comparison of our radiative transfer 
modeling with these observations. This is necessarily a crude 
comparison because we are limited to a single quasar (albeit 
simulated at several luminosities) situated in a dark matter 
halo which is most likely less massive than those hosting 
the quasars in the HP06 sample. In Figure \T2[ we show the 
distribution in luminosities for the Nios =61 foreground 
quasars in their sample whose sight lines have impact pa- 
rameters that would place them inside our cut out proximity 



iHennawi et"al] ([2006^ report the maximum enhance- 
ment of the UV background due to the foreground quasars 
in their sample as, 

Fqso 



guv = 1 + 



(10) 



where Fqso is the number flux of ionizing photons at the 
position along the line of sight closest to the foreground 
quasar and Fjjv is the number flux of photons at the same 
position due to the UV background. To translate this into 
an upper limit for the number of ionizing photons emitted 
by the foreground quasar per second we use. 



Nqso — 'inb^ Fqso — 47r6^ Fuv{guv 



(11) 



where b is the impact parameter of the line of sight to the 
background quasar. We have used a spectral slope of q = 1.8 
and a luminosity density at th e Hydrogen ionizing thre sh- 
old of Jhi = 5.0 X 10"^^ from iHaardt fc Madaul |l99i) to 
calculate Fuv- 

Also shown in Figure [12] is the mean value of the HP06 
sample luminosities (black horizontal solid line) and the 
seven quasar luminosities we used in our simulations (blue 
horizontal dashed lines). The vertical cyan lines in the plot 
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Table 2. Neutral Hydrogen Statistics. 



Radiative 
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Shielding 
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Model 
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NA 
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Figure 11. MoUweide projections showing the fluctuation of Log Ffji on spherical shells at several radii. From left to right and top to 
bottom, the radii are R = 12.5, 25, 50, 100, 200, and 400 kpc/h from the central black hole. These are the same radii that are shown in 
Figure [To] and the fluctuations are around the peak at each radius in that plot. A persistent shadow can be seen in the lower left region 
of each plot. 
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indicate the impact parameter bins we used to construct a 
sample of A'mock = 1000 sets of Nios ~ 61 lines of sight with 
the same distribution in impact parameter as their sample. 

In Figure [131 we present our comparison. For the fidu- 
cial ray tracing model, we count the number of lines of 
sight above threshold in each of the A'mock groups and his- 
togram the values. We repeat this procedure for each of the 
seven simulated quasar luminosities (see Figure I12|) and plot 
the results as thin yellow lines in the upper left (Log 
= 17.2) and lower left (Log N^^ =19.0) panels. Overlaid on 
these lines in thick black is the luminosity weighted aver- 
age of each bin. The weighting is done by determining the 
fraction of HP06 quasars that is nearest to each simulated 
quasar luminosity. In the second and third columns we re- 
peat this procedure for two semi analytic models. The last 
two columns show the same thing for ray tracing runs that 
include only a UV background sent in from outside the box 
(fiducial in fourth column and half the fiducial in fifth col- 
umn). These two background intensities produce photoion- 
ization rates within ~ 3(7 of the value determined by FG08 
(see Figure ID). We will use these panels to discuss the effects 
of obscuration and periodic emission shortly. We have indi- 
cated the observed values from HP06 with vertical dashed 
lines and the peak of the luminosity weighted average from 
the first panel ( GADGET -I- SPHRAY ray tracing runs) with 
arrows in each panel. In all cases which include contributions 
from the central quasar, we fall short of the observed inci- 
dence of optically thick absorbers. There are several plausi- 
ble explanations for this. 

Before we present these explanations we briefly discuss 
the pixel maps used to generate Figure 1131 First, we note 
that the semi analytic model that comes closest to repro- 
ducing the observational results is unrealistic in terms of 
the distribution of HI. In Figure [141 we show the projected 
column densities and the positions of galaxies identifled in 
the base hydrodynamic simulation. The top row is the di- 
rect GADGET output, the top middle row is the fiducial ray 
tracing GADGET + SPHRAY run, the bottom middle row 
is a semi analytic model with no shielding and the bottom 
row is the Log A = 4.0 semi analytic run. The columns show 
magnifications of the central region. 

7.1 Proximate Density Field - Host Halo 

The most straightforward explanation for the disagreement 
we see is that our single density field is not representative 
of the proximity zones in the IIP06 sample. The neutral 
Hydrogen density field is determined by a combination of 
the total baryonic overdensity in the proximity zone, the 
background radiation field, and the ionizing flux from the 
central point source. We have taken account of the various 
quasar luminosities by running multiple models, however we 
are restricted, in this study, to a single density fleld. 

The mass of the dark matter halo hosting our central 
quasar is MhsIo = 5.25 x 10^^ h~^M(^. Observations of the 
clustering of qua,sars (e.g. Porciani et al.l 20041: Groom et al.l 



l2005l : IShen et al.ll2007l : ICoil et al.ll2007l : rShen et al.l2009D im- 
ply a minimum quasar host halo mass between lO^'^ and 
lO" /i~^M0 with a typical mass closer to lO" /i'^Mq. In 
the third of the Quasars Pro bing Quasars series of papers, 
jProchaska fc Hennawil (|2009l ) examine a single representa- 
tive quasar pair in detail. In this analysis, they use a mass 



of A/Halo = 2.01 X 10^^ /i~^Mq for the foreground quasar's 
host halo mass. Therefore, it is plausible that our dark mat- 
ter host halo is a factor of four less massive than the typical 
host halo from th e HP0 6 sample. 

iKim fc CroftI | 20081 ) have constrained host halo masses 
at a mean redshift of z = 3 using the mean transmitted 
flux in the Ly-a forest and find agreement with the previ- 
ous works mentioned above. In addition, they have shown 
that the baryonic over density within a given radius is an 
increasing function of halo mass. Specifically, in Figure 2, 
they show an increase in baryonic overdensity by a fac- 
tor of 3.66 as dark matter halos range from 1.6 x 10^^ to 
5.1 X 10^^ MqH-^ using a 100 h '^Mpc cosmological sim- 
ulation. Our fiducial ray tracing run peaks at a value of 5 
absorbers (top left panel in Figure [TSJ with column densities 
Log A^Hi ^ 17.2. HP06 observe 15 lines of sight with impact 
parameters that would put them in our cut out region above 
this threshold . Usin g a simple linear scaling from the data of 
iKim fc CroftI (l2008h in the relevant mass range and for the 
number of absorbers with overdensity, we estimate the rep- 
resentative dark matter host halo mass in the HP06 sample 
as M = 1.92 X 10^^ Mf^h''^ . This agr ees well with the value 
used in lProchaska fc Hennawil l|2009l ). 



7.2 Unresolved Clumping 

Aside from the fact that our theoretical sample is limited 
to a single density field, the dynamic range required to re- 
solve small scale clumping of gas in a cosmological volume 
is extreme. In the base hydrodynamic simulation, the dis- 
tribution of smoothing lengths peaks at « 8.5 h~^ kpc over 
the whole simulation, and ~ 850 h~^ pc or 10"* of the full 
box size in dense regi o ns. In a follow up paper to HP06, 
IProchaska fc Hennawil l|2009l ) estimate the size of a single 
absorber as 10-100 pc while the smallest smoothing length 
in our cutout proximity zone is 150 h^^ pc. Better resolv- 
ing this small scale clumping would increase the incidence of 
optically thick absorbers and increase the column density of 
those already identified. While it would be a very large un- 
dertaking to simulate an entire cosmological volume at the 
necessary resolution, it is feasible to re-simulate a handful 
of proximity zones at higher resolutions. On the other hand, 
we could enlarge our cut out region further and use semi- 
analytic treatments of this clumping by requiring our global 
(as opposed to that around a single quasar) distribution of 
Damped Lyman-a systems to match observations. In prac- 
tice, simulations like these are used together to calibrate one 
another and both are something we plan to do in the future. 



7.3 Obscuration and Duty Cycle 

In unified theo retical models of AGN and quasars (e.g. 
lAntonuccilll993l '). the emission from the accretion disc and 
broad line regions is obscured over a certain solid angle by 
an optically thick torus. The inhomogeneities this could in- 
duce in the Lyman a forest have b een theoretically inves- 
tigate d bv iGroft (2004 ). In addition. iHennawi fc Prochaskal 
((2003) find evidence of this obscuration by comparing the 
incidence of optically thick absorbers in foreground quasar 
lines of sight with the transverse lines of sight to background 
quasars. Our simulated sources emit isotropically and so. 
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even with self shielding, there could be large solid angles in 
which we are overestimating the photoionization rate. To get 
an idea of the magnitude of this effect we present the back- 
ground only panels in Figure [131 (columns 4 and 5). These 
backgrounds bracket the value calculated in FG08 and are 
both within 3 a of the backgrounds present at the bulk of the 
redshifts in the HP06 sample. If the radiation is only present 
within certain solid angles, the incidence of absorbers in our 
single density field will lie somewhere between the ray trac- 
ing panel and one of the background only panels. This could 
account for a non negligible part of the discrepancy we see 
as 1.8% and 5.3% (columns 4 and 5) of the mock groups of 
sight lines overlap the HP06 data. 

In a similar way, a central black hole undergoing an 
episode of accretion that lasts less that the light crossing 
time of the region considered can lead to large volumes 
where we have overestimated the photoionization rate. The 
light crossing time for our cutout region, from the center 
of the cube to the corner, is tx ~ 6 Myr. The lower limit 
on episodic emission from quasars is set by observations of 
the line of s ight proximity ofi'cct in the Ly-a forest (e.g. 
iBaitlik et all ,1988: Scott ct al. 200(1) and is equal to the 
photo ionization equilibration tim e of the IGM, teq = 10* 
years l|Martini fc Weinberg|[200ll ). If the central source were 
on for a time ton where 0.01 ^ ton ^ 6 Myr then certain 
radii would not have felt the enhanced UV field from the 
central source yet. 



7.4 Velocity Offsets 

It is straightforward to determine the angular separation 
of two objects on the sky. Due to redshift distortions, it 
is much more difficult to determine their separation in the 
radial direction. This problem is exacerbated for gas near 
large overdensities where it has large infall / outflow veloc- 
ities. The definition of proximate DLA in the HP06 sample 
is systems with absorber redshifts within 3000 km/s of the 
emission redshift of the foreground quasar. Due to this un- 
certain radial separation, we have included all absorbers in 
the observational sample that satisfy the impact parameter 
requirements of our cut out region. This includes some ab- 
sorbers from the HP06 sample in our comparison that would 
not actually be in our cutout region. A larger cut out region 
would allow us to better model these errors. 



8 CONCLUSION AND DISCUSSION 

In this paper, we have presented a detailed model for the 
neutral hydrogen gas in the proximity zone of a quasar pow- 
ered by accretion onto a super massive black hole. To do 
this, we combined a cosmological hydrodynamic GADGET 
simulation of the formation and growth of black holes with 
a detailed radiative transfer post processing code SPHRAY. 
We focused our attention on self shielding of the neutral 
gas and, by construction, the contribution from a local cen- 
tral source. Much of the previous work on Lyman Limit and 
Damped Lyman a systems ignored the contribution from 
point sources assuming that the background field dominates. 
This is true, on average, for absorbers less dense then Lyman 
Limit systems but is not true for the more dense absorbers 
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Figure 12. Maximum foreground quasar luminosity vs impact 
parameter for the sample of quasar pairs in HP06. Pairs that 
contain an absorber with column density LogAfni > 17.2 are col- 
ored red while pairs that contain an absorber with column density 
LogTVni > 19.0 are colored green. We have indicated the mean 
of the observed luminosities with a black line and the luminosi- 
ties of the central black hole from our simulations in dashed blue 
lines. The vertical cyan lines indicate the bins used to calculate 
a set of mock lines of sight with the same distribution in impact 
parameter as the observations. 



which cluster around luminous point sources and those un- 
derdense systems that just happen to be near sources. In ad- 
dition it is crucial that these sorts of models be constructed 
to interpret quasar pair observations such as those in HP06. 

We find that the results of our ray tracing algorithm 
SPHRAY cannot be reproduced with simple semi analytic 
models in which gas above a given density threshold is com- 
pletely shielded from ionizing radiation. This is because the 
photoionization rate at any given point is determined by the 
partial shielding provided by the local density field and to a 
lesser degree the partial shadowing seen in Figure [TT] 

We have shown in Figure |9] that our GADGET -|- 
SPHRAY models naturally reproduce the observed shape 
of the A'hi column density distribution. We have also 
shown that the flattening feature around Log A''hi = 
20 attributed to self sh i elding in the analytic work of 
I Zheng fc Miralda-Escudi l|2002l i is not due entirely to self 
shielding, but is largely enhanced by it. 

When compared to the observations of HP06, our ray 
tracing and semi analytic models fall short of capturing the 
number of dense absorbers proximate to luminous quasars. 
We have described improvements to our model that would 
bring our results into agreement with these observations in- 
cluding sampling more massive host halos, correcting for the 
under resolved clumping of gas, allowing obscuration / pe- 
riodic emission of the quasar radiation and using a larger 
simulation volume to better model uncertain radial separa- 
tions. These issues are not intractable and can be dealt with 
using a larger low resolution simulation to sample more mas- 
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Figure 13. Lines of siglit with column densities A'^jjl above for several radiative models. We construct Nj^^^]^^ = 1000 groups of 

^laa = 61 lines of sight that match the impact parameter distribution of the HP06 sample. For the fiducial ray tracing model, we count 
the number of lines of sight above threshold in each of the Af„jQ(.]j groups and histogram the values. We repeat this procedure for each of 
the seven simulated quasar luminosities (see Figure Fl2l l and plot the results as thin yellow lines in the upper left (Log N^^ =17.2) and 
lower left (Log =19.0) panels. Overlaid on these lines in thick black is the luminosity weighted average of each bin. The weighting is 
done by determining the fraction of HP06 quasars that is nearest to each simulated quasar luminosity. In the second and third columns 
we repeat this procedure for two semi analytic models. The last two columns show the same thing for ray tracing runs that include 
only a UV background sent in from outside the box (fiducial in fourth column and half the fiducial in fifth column). These are useful 
in describing the effects of beaming. We have indicated the observed values from HP06 ( vertical dashed lines) and the peaJc of the 
luminosity weighted average from the fiducial ray tracing runs (arrow) in each panel. 



sive dark matter host halos coupled with a high resolution 
resimulation of a more representative proximity zone. 

Finally we have shown that the total amount of neu- 
tral HI within the 750 kpc region around our central 
quasar can change by up to an order of magnitude depend- 
ing on the shielding prescription used. This means that HI 
surveys that do not resolve individual galaxy features but 
whose goal is to measure the integrated signal from a prox- 
imity zone will need to take these effects into account. This 
is on sca les much smaller then the intensity mapping pro- 
posed bv lPeterson et all (|2009l ) to measure Baryon Acoustic 
Oscillations, but would be relevant for connecting HI or 21 
cm galaxy surveys to the underlying density field. 
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APPENDIX A: RADIATIVE TRANSFER 
CONVERGENCE 

Here we examine the convergence of our ray tracing radiative 
transfer runs. In these runs, both the photoionization rate 
Phi and the neutral fractions xm are taken as outputs from 
SPHRAY . In Figure [All we show the evolution of xm as the 
number of rays traced is increased over four decades. The 
distribution has a triple peaked structure with the left peak 
being due mostly to collisional ionizations, the central peak 
due mostly to photoionizations in the diffuse IGM and the 
right peak formed by self shielded particles. The photoion- 
ization peak has not converged completely however there is 
a very small amount of neutral material here and our main 
interest is in the self shielded peak which has very nearly 
reached convergence. 

In Figure IA2I we examine the convergence of Fhi • For 
studying self shielding, the position of the peak is not as 
important as the tail of the distribution extending to smaller 
photoionization rates. The slightly smaller values in the tail 
for the rt9 run translate into slightly less neutral material 
with —3 < Xm < — 1- 



REFERENCES 

Altay G., Croft R. A. C, Pelupessy I., 2008, MNRAS, 386, 
1931 

Antonucci R., 1993, ARA&A, 31, 473 
Bagla J. S., 2002, JApA, 23, 185 

Bajtlik S., Duncan R. C, Ostriker J. P., 1988, ApJ, 327, 
570 



14 Gabriel Altay et al. 



IB 



23 



ai 17 



1B 



10 



l^ 17 



1H 



aa 



ai 



n- : ■ 


■ T ,1-1 #j t * ■ 












alt. 






















* 



























Si 







-BDQ-^iO iC3 3 JGD KB BDD 



1D!X D 



-a 3 ZD 

K EWi] 



Figure 14. The HI column density obtained by collapsing the 1500 kpc/h simulation volume along the z-axis. The top row shows the 
direct GADGET output with no central UV field from the quasar, the top middle row shows the GADGET + SPHRAY model, the 
bottom middle row shows the semi analytic model with no shielding, and the bottom row shows a semi analytic model with a shielding 
threshold at A = 4.0. The central quasar has the same luminosity in all three models which include it. Overlaid in solid cyan lines are 
the galaxies that were identified in the parent hydrodynamic simulation. The radius of each circle is proportional to the cube root of the 
gas mass. We have indicated Log M[M0/h] for three galaxies chosen to be representative of the spread in mass. The symbol with Log 
M = 11.59 is the central quasar host galaxy. 
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Figure Al. Mass weighted probability distribution of xhi for 
the fiducial GADGET + SPHRAY ray tracing model including 
the background and the central quasar with varying numbers of 
rays traced. Shown as an inset is a close up of the self shielding 
neutral peak. 
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Figure A2. Mass weighted probability distribution of Fhi for 
the fiducial GADGET + SPHRAY ray tracing model including 
the background and the central quasar with varying numbers of 
rays traced. Shown as an inset is a close up of the self shielding 
neutral peak. 
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